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^3 ' Abstract 

Age structure is incorporated in many types of epidemic model. Often it is convenient to 
assume that such models converge to early asymptotic behaviour quickly, before the susceptible 
00 ' population has been appreciably depleted. We make use of dynamical systems theory to show 

that for some reasonable parameter values, this convergence can be slow. Such a possibility 
should therefore be considered when parameterising age-structured epidemic models. 



1 Introduction 



d 

I 1 Age is one of the key variables to consider in any epidemiological analysis, particularly related 

, ^ , ', to infectious disease (Anderson and May, 1991). In the context of human diseases, age is an 

extrinsic risk factor, strongly influencing an individual's mixing patterns within the population 
(Keeling and Rohani, 2008); but age is also an intrinsic risk factor for serious outcomes and this 
can lead to important and often counter-intuitive public-health conclusions (Anderson and May, 
1983). 

While classical models of childhood infections used a priori matrices to parameterise the epidemio- 



logical mixing between individuals of different ages (Anderson and May, 1991; Keeling and Rohani, 



2008), more recent studies have attempted to measure this quantity empirically either using surveys 



(Edmunds et al., 2006; Mossong et al., 2008) or through the use of synthetic populations derived 
from large datasets (del Valfe et al., 2007). The insights gained by doing this were particularly 
useful during the 2009 influenza pandemic (Baguelin et al., 2010). A mathematical overview of the 
impact of recent empirical results for modelling is given by Glasser et al. (2012). 



In theoretical analysis of age-structured models, it is often convenient to make the assumption that 
early behaviour of the system has converged to an appropriately defined dominant eigenvector, 
representing the relative prevalence of disease in different age groups. Often this is done through 
the consideration of discrete generations of infectious individuals, and the relevant eigenvector is 
associated with a 'next-generation matrix' whose dominant eigenvalue is the basic rcoproductive 
ratio Rq, which is the expected number of secondary cases per primary case in a naive population 
(Dickmann et al., 1990). While discrete generations are not always straightforwardly related to the 
system's real-time behaviour, this framework opens up the possibility of highly general analysis of 
a large number of epidemiological scenarios (Diekmann and Heesterbeek, 2000). 

Here, we analyse the rate of convergence of an age-structured epidemic to its dominant eigenvalue in 
a real-time framework, using both analytic and numerical methods. This is done using a dynamical 
system of ordinary differential equations (ODEs), together with results from linear algebra. We 
find that for some plausible parameter values, it is possible for the timescale of convergence to 
be comparable to the timescale over which susceptible depletion becomes appreciable, and so the 
assumption of fast convergence to the dominant eigenvector cannot be made. 
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2 Age-structured model 



2.1 Model definition 

Our modelling approach is based on the SIR compartmental structure in which individuals are all 
cither Susceptible, Infectious or Recovered, without demographic processes like birth and death. 
Individuals arc also placed into discrete age categories indexed by a, & € N, rather than using con- 
tinuous ages that would require a less tractable partial differential equation or intcgro-diffcrential 
equation model. Our ODE system is then 

-jjT = Sg ^ Pa.bh - jla ■ 
b 

Here the transmission matrix is M = (/3 ,b)j where (3 a ^ is the rate of transmission from individuals 
in age class b to individuals in age class a. 

We will assume that there are n age classes and adopt, for convenience, a normalisation 

S a (t) + I a (t) + R a (t) = 1 , Va,t. (2.2) 

This simplifies the algebra involved in analytic work (i.e. a more general case can be considered in 
our framework at the cost of more complex and less enlightening expressions) and can be justified 
for particular choices of age classes. 

Primarily our approach will look at analytically approximating the solutions of our system near 
fixed points to determine the behaviour of infecteds and susceptibles in the population. We will 
adopt a vector notation to represent the dynamical state of the system. 

/ Si \ I h \ 

In this notation, we consider dynamical systems of the general form 

f=F(x). (2.4) 

If x* is a vector such that F(x*) = 0, the system dynamics can be linearised around this vector to 
give 

dy 

— = Jy + 0(e) , where x = x* + ey . (2.5) 

In this regime, the dynamics are dominated by the eigenvalue of the Jacobian matrix J with the 
largest real component. We call this the dominant eigenvalue, written Ai, which has an associated 
dominant eigenvector, v±. Typically, there will be a rate r such that the ratio of the magnitude 
of y in the direction of vi to the magnitude of y in any orthogonal direction will asymptotically 
grow like e rt . We will have a quantity like r in mind when discussing the rate of convergence of the 
system to its dominant eigenvalue. If this rate is low compared to the rate at which 0(e) effects 
become important, then we will say that the system fails to converge to its dominant eigenvector. 
We note that this is different from the rate of convergence of a stochastic model to its deterministic 
limit (Tuljapurkar, 1982), and is fundamentally a real-time concept, meaning that direct analogues 
need not exist in a discrete generation-based framework as analysed by Dickmann and Heesterbcck 
(2000). 
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2.2 Early dynamical behaviour 



In order to consider behaviour early in the epidemic, the Jacobian of (2.1) will have to be considered, 
and is given in terms of quantities defined in the previous section by 
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We are interested in the behaviour around the disease-free equilibrium of (2.1), x* = (1, 0) T , where 
1 and are a length-n vectors whose entries are all 1 and all respectively. Putting this fixed 
point into the Jacobian gives 



J(x*) 



-M 
M-G 
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(2.6) 



Then the early behaviour of the epidemic is given by setting 1(0) = + ey, leading to early 
dynamics of the form 

±I(t) = (M-G)I(t) + 0(e 2 ) . (2.7) 



2.3 Perron-Frobenius analysis 

In general, we expect the matrix M — G to be quasi-positive, since we expect every age to interact 
with every other at some level, leading to positive off-diagonal elements, but recovery can lead to 
elements < on the diagonal. We consider the special case where each age class is capable of 
supporting the disease independently of the others (i.e. fi a . a > 7) which is more likely to apply 
if the disease is more transmissible and there are fewer age classes in the model. In this special 
case, the matrix M— G has positive real entries, allowing us to use the following important theorem. 

Theorem 1. (Perron-Frobenius) Let A = (ajj) be an nxn positive matrix: a^j > for 1 < i,j < 
n. Then the following statements hold: 

i. There is a positive real number Ai, called the Perron root or the Perron-Frobenius eigen- 
value, such that Ai is an eigenvalue of A and any other eigenvalue Ai^i has \Xi\ < X\. 

ii. The Perron-Frobenius eigenvalue is simple. 

in. There exists an eigenvector vi of A with eigenvalue X\ such that all components of Vi are 
positive. 

iv. There are no other positive eigenvectors. I.e. all other eigenvectors must have at least one 
negative or non-real component. 

Proof. See e.g. Meyer (2000, Chapter 8). □ 
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We will also assume that the initial vector direction for the dynamical system, y, can be written 
as a linear combination of the eigenvectors {v,} of M — G 

n 

y = ^C ?; v,, (CieC). (2.8) 
This, together with (2.7), means that early in the epidemic, 

n 

1(f) =£c i v i e Ait , (ceC) 

i=l 

\ i=2 

From the Perron- Frobenius theorem, we can then argue that for sufficiently small |cj>i|, and 
sufficiently large t, the contribution to I from eigenvectors other than vi can be made arbitrarily 
small since the absolute value of such contributions decays exponentially over time. 




3 Analytical approach 

3.1 Models of assortative mixing 

We now consider mixing matrices that have the assortative property observed in real networks, but 
remain amenable to analytic methods. These are in fact special cases of the matrices considered us- 
ing a time- independent framework to determine invasion thresholds in Diekmann and Heesterbeek 
(2000, §5.3.2). Our aim, however, is to consider time-dependent transient behaviour for diseases 
that are successfully invading. 



Definition 2. A Basic Mixing Matrix is an n x n matrix B 



B(a, n) 



( 1 + Q 
1 

1 



V i 



f 

1 + Q 
1 



1 



a 



such that 

1 
1 



1 

1 + 



+ a J 



Where a G R, a > — 1 . 

Theorem 3. Let B(a, n) be a Basic Mixing Matrix, then A = n+a is an eigenvalue with eigenvec- 
tor vi = (1,1,..., 1) T and A = a is an eigenvalue with algebraic multiplicity n — 1 corresponding 
to n — 1 eigenvectors of the form; 



Vj^i = (1,1, -n,...,l) 



with 1 — n in the i th place. 



Proof. Eigenvalues and eigenvectors satisfy the equation Bv = Av, hence letting vi = (1, 1, . . . , 1) T 
gives A = Ai = n + a, and by Perron-Frobenius, this eigenvalue is simple and all other eigenvalues 
are less than A with eigenvectors that have at least one negative component. If we now take 
Vj = (1, 1, . . . , 1 — n, . . . , 1) T , i = 1, . . . , n for i ^ 1, where 1 — n is in the i th place, this gives 
A, = a, Vi 7^ 1. Then since Vj are linearly independent Vi we have found the entire system of n 
eigenvectors and A = A,^i = a has algebraic multiplicity n — 1. □ 
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Definition 4. A Simple Mixing Matrix is an n x n matrix S such that 



( 1 + ai 1 
1 l + a 2 



S(a,n) 



1 



V i 



1 



1 



1 
1 

1 + a 3 



1 1+ttn / 



Miere a t G R, a, > — 1. 

Theorem 5. TTie Simple Mixing Matrix S(a,2) ftos eigensystem 



2 + ai + a 2 ± v(ai - 02) 



v b = 



CK2— aii-^/ (ai — 02) 2 +4 
2 



Proof. Substitute (3.1) into S(a,2)v± = A±v±. 



(3.1) 



□ 



Note that finding the eigensystem of S(a, n) in general involves solving an order-n polynomial, but 
no simple general formula exists as for a Basic Mixing Matrix, limiting the usefulness of analytic 
methods, and we introduce this concept mainly to explain why a more general analytical approach 
is not possible. 



3.2 Early time estimate 

We now wish to define a timescale over which the epidemic approaches its dominant eigenvalue. 

Definition 6. Consider an epidemic model where the early vector of infectious individuals is given 
by (2.9). The Early Time is defined as 



1 



A- A 



where A := max {Rc(A;>i)} 



(3.2) 



This definition has the property that over a period of, say, time 2r, the relative importance of each 
non-dominant eigenvector reduces by a factor of at least e 2 w 7.389. 



3.3 Late time estimate 

An estimation of the late time is also necessary to assess if convergence to the dominant eigen- 
value happens before susceptibles are depleted. Making use of the standard definition of Rq 
(Dickmann et al., 1990) as the dominant eigenvalue of the next-generation matrix, in our case 
M/7, we expect intuitively that the epidemic is near to its maximum when the proportion of the 
population susceptible is S(t*) = I/Rq, since for this susceptible population (if allocated appropri- 
ately to different age classes) the disease can no longer invade. 

Definition 7. Consider an epidemic model with basic reproductive ratio Rq given by the dominant 
eigenvalue of M/7, starting with a proportion s of the population infectious. The Late Time is 
defined as 

t -=H^)- (3 - 3) 

This definition makes sense since if the population is initially entirely susceptible we may ap- 
proximate S(t) w 1 — ec A *. The time at which this quantity first exceeds l/i?o should therefore 
correspond to a sensible definition of late time. Figure 1 visualises the concepts involved in this 
definition. 
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3.4 Separation of early and late time 



We are now in a position to compare early and late time. Clearly, it is always possible to make e 
sufficiently small to keep these timcscales separated; explicitly, this requires 

T<<T ^ e «i^ exp ( I A_). ( , 4) 

Now suppose we consider transmission based on a Basic Mixing Matrix multiplied by an overall 
transmission rate v so that M = i/B(a, n), meaning that the eigensystem of M — G can be obtained 
by application of Theorem 3. The basic reproductive ratio is Rq = v(n + a)/7, and for an initial 
infection of Io in the i th age group, 

I(0) = (0,...,0,/o,0,...,0) T = ^ (vi-Vi) . (3.5) 
i— i n 

i-th position 

Note that the symmetry of the system means that we can let i ^ 1 without loss of generality. This 
initial condition leads to early behaviour 

10) = ^c At (vi - e~ (A_A)t v^ , where A = v{n + a) - 7 , A = z/a - 7 , (3.6) 

and so the early time r = l/(yn) . It is most instructive to consider what happens at constant Rq, 
so substituting into (3.4) in this regime gives separation of timescales when 

Io < —5 — exp — — . (3.7) 

R v \y a h) - R oJ 

Taking appropriate partial derivatives of the right-hand side of this expression, we see that for a 
given Rq, a lower Iq is required to secure convergence if we increase 1/ or a. or if we decrease 7. 
This is as would be expected from thinking qualitatively about the problem. 



4 Numerical approach 



4.1 Quantification of convergence 



It will be convenient to define a numerical measure <f>{t) of convergence to the dominant eigen- 
value. 

4>{t) := tan" 1 ^)) , 



d / net) 

where ip(t) := - In' 



dt \ \J{t) 
I(t) :=J2 T S) , 



(4.1) 



J(t) := e At . 

This definition gives <p{t) = at times when the epidemic moves along the dominant eigenvector, 
and non-zero values in the interval (— 1, 1) otherwise. 



4.2 Contact survey data 

We take the form for M in our numerical work from the results of the POLYMOD study (Mossong et al., 
2008) for all reported contacts (physical and conversational) across all countries. We use age 
classes of width 5 years and a final class for ages 70+. Consideration of population pyramids 
(Wroth and Wiles, 2009; Stillwell and Clarke, 2011) shows that each age group is of almost equal 
weighting up to 70, after which the rest of the population is in a rapidly declining tail which in 
total is roughly equal to the weighting of one of the 5 year classes, and so for present purposes 
where we are looking for broad insights rather than precise predictions the normalisation (2.2) is 
appropriate. 



G 



4.3 Results 



When we integrate (2.1) directly using Runge-Kutta, for infection starting in the 15-19 age group, 
Figure 2 shows the results obtained for cf>(t) and I(t). Strikingly, this plot shows that even during 
periods of the epidemic where the prevalence timeseries looks linear on a logarithmic y-axis, 4>(t) 
can fail to equal zero for an appreciable period of time (meaning that the dynamics (2.1) are not 
governed by the dominant eigenvector) for modest initial infectious populations. Whether such 
initial concentrations of infection can be seen, perhaps as a result of an unpredictable event early 
in the epidemic when stochastic effects dominate, depends on the population size that can be seen 
as mixing according to (2.1). At country size - i.e. tens or hundreds of millions of individuals - 
small Io is most likely, but at city sizes of around 10 5 , it is not inconceivable that an epidemic 
could start with several hundred students in the same school leading to Io ~ 10~ 3 . 

5 Discussion 

In this paper, we have considered the real-time rate of convergence of age-structured SIR epidemic 
models to their dominant eigenvector, both analytically for a basic model of assortative mixing, 
and numerically for realistic mixing. In each case, sufficiently small initially infectious populations 
arc needed to ensure that the dynamical system does indeed converge before the depletion of 
susceptibles becomes important. In many cases this effect will not matter; but it is plausible 
that for some supercritical epidemics in moderately-sized, closed populations with strongly age- 
structured mixing, this failure of convergence to the dominant eigenvector will be a source of 
bias. In particular, attempts to infer the relative susceptibility of different age classes may become 
biased. If a particular age group is over-represented in the initially infected group as compared 
to the dominant eigenvector, then its relative susceptibility may be overestimated. Correction for 
non-convergence to the dominant eigenvector is therefore potentially important in epidemiological 
analysis. 
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Figure 1: Visualisation of the late time approximation. Peak time t* corresponds to the peak of 
the full epidemic (solid black line) and late time T is the equivalent timescale for the exponential 
growth approximation (dashed line). 





Figure 2: Top: Distance from the dominant eigenvalue, 0, for the POLYMOD mixing matrix. 
Bottom: prevalence against time on a logarithmic scale. Parameters arc: i?o = 5, 7 = 0.1, 
I Q = lCT 1 (solid lines), I = 1CT 3 (dashed lines), J = 1(T 5 (dash-dot lines), I = 1CT 7 (dotted 
lines). Infection is started in the 15-19 year-old age group. 
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